A quantum effect in the classical limit: 
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For suitable parameters, the classical Duffing oscillator has a well-understood bistability in its 
stationary states, with low- and high-amplitude branches connected by an unstable intermediate 
branch. As expected from the analogy with a particle in a double-well potential, transitions be- 
tween these states become possible either at finite temperature due to noise, or in the quantum 
regime due to tunneling. In this analogy, besides local stability to small perturbations, one can 
also discuss global stability by comparing the two potential minima. For the Duffing oscillator, the 
stationary states emerge dynamically from an interplay between driving and dissipation so that a 
priori, a potential-minimum criterion for them does not exist locally, let alone globally. However, 
global stability is still relevant, and definable as the state containing the majority (in the limit, all) 
population for long times, low temperature, and close to the classical limit. Further, the crossover 
point is the parameter value at which global stability abruptly changes from one state to the other. 

For the double-well model, the crossover point is unambiguously defined by potential-minimum 
degeneracy. Given that this analogy is so effective in other respects, it is thus striking that for 
the Duffing oscillator, the crossover point turns out to be non-unique. Rather, none of the three 
aforementioned limits commute with each other, and the limiting behaviour depends on the order 
in which they are taken. More generally, as both h — > and T — > 0, the ratio huio/k&T continues 
to be a key parameter and can have any nonnegative value. This points to an apparent conceptual 
difference between equilibrium and nonequilibrium tunneling. We present numerical evidence for 
this phenomenon by studying the pertinent quantum master equation in the photon-number ba- 
sis. Independent verification and some further understanding are obtained using a semi-analytical 
approach in the coherent-state representation. 



PACS numbers: 03.65.Yz, 05.40.-a, 85.25. Cp 

Introduction. The Duffing model describes an oscilla- 
tor with weak damping, driving, and a quartic potential 
nonlinearity lj. Given that a cubic nonlinearity renders 
the potential unbounded from below on one side, with 
essentially different physics, and that the quartic term is 
the next one encountered in the expansion of a generic 
one-dimensional potential, the model is a prototype for 
weakly anharmonic oscillators. As such, it is unsurpris- 
ing that it is encountered when describing a wide range 
of oscillation phenomena, be they mechanical, electric, 
or optical. For a range of parameters (discussed in detail 
below), the model exhibits a well-known classical bista- 
bility. This leads to the application serving as our main 
motivation, as a superconducting readout device for su- 
perconducting quantum bits (qubits) [2]. In this context, 
the Duffing oscillator is known as the Josephson bifurca- 
tion amplifier (JBA) 3 . A key advantage of the JBA is 
that both bistable stationary states are superconducting, 
so that one avoids the power dissipation and quasiparti- 
cle generation which would accompany a switch to the 
normally conducting state. However, our focus will be 
on the oscillator proper, not on its coupling to a qubit or 
other device. 

In the JBA application, the Duffing oscillator is suffi- 



ciently large to be in the classical regime, and a classi- 
cal description suffices for the main experimental results. 
Nonetheless, a fundamentally consistent treatment of the 
coupling to a quantum device necessitates a quantum- 
mechanical description of the entire system [3]. Further, 
states which are (bi)stable classically may acquire a fi- 
nite lifetime due to quantum tunneling, which should be 
properly quantified even if it is negligible in a given exper- 
iment [5HZ|- More ambitiously, one could try to design 
an experiment focusing on this unconventional form of 
nonequilibrium tunneling in a driven system, which will 
indeed turn out to have some unusual and novel features. 

First, we will review the classical situation. Subse- 
quently, a quantum formulation is developed and an- 
alyzed numerically, introducing and characterizing the 
phenomenon of crossover. As in other contexts, quantum 
tunneling in the Duffing oscillator can be contrasted with 
thermal activation and diffusion at finite temperature, 
which persists in, though is not restricted to, the clas- 
sical limit. Comparison of the two approaches reveals a 
striking difference in crossover behaviour in the classical, 
low-temperature limit. Therefore, we focus on this limit 
using semi-analytical asymptotic methods. In particu- 
lar, it is physically transparent (though not convenient 
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numerically) to cast the quantum equations in a form 
closely resembling the classical ones, using a coherent- 
state representation. We finish by making some conclud- 
ing remarks. 

Classical dynamics. For background, and to establish 
some notation, let us recapitulate the classical situation. 
The Newton-Duffing equation is 



H' = 3.24 



x + uj 2 x = — e[yx + ax 3 + f cos(wi)] 



(1) 



in which the driving strength /, the damping 7, and the 
potential nonlinearity a [all defined such that the mass m 
cancels in Eq. |l])] are governed by the same small param- 
eter e —> 0. Due to the oscillator's high quality, nonper- 
turbative physics can nonetheless emerge for long times 
if the driving is sufficiently close to resonance, 



(2) 



with Q denoting the detuning. Since the potential an- 
harmonicity is ^emax 4 , we take a > to have a globally 
stable potential |S]. 

To focus on the long-time dynamics, transform the fast 
(~ Uq 1 ): large, near-harmonic oscillations away by going 
to the rotating frame, in this context also known as Van 
der Pol coordinates [9]: 



u\ I cos uit — sinwA / x 
v! \ sin ujt coswy \p/uim 



(3) 



with p = mx. One finds that the expressions for (u, v) in- 
deed are 0(e), so that (it, v) are varying slowly; thus, the 
integrated effect of terms such as wcos(wi) in the equa- 
tions of motion is of higher order in e, hence negligible. 
After this coarse-graining, the surviving slow dynamics is 
expressed in its natural (dissipative) time scale r = eft, 
as 



d_ fu 
dr \V 



Ju 
Jv 



with the scaled rotating-frame flow field (see Fig. [IJ 
1 ( -fi'V -U + \a'V{U 2 + V 2 ) \ 



J = 



il'U -V-l- \a'U(U 2 + V 2 ). 



and the dimensionless variables If 2 



Q' 



w 7 



af 2 



^o7 
/ 



(4) 



(5) 



(6) 



In particular, detuning is scaled relative to the natural 
linewidth, while coordinates are relative to the resonant 
amplitude. 

Unlike the driven Eq. ([!]), which is equivalent to an 
autonomous three-dimensional system and is known to 
have a chaotic regime jTU] , Eq. Q is an autonomous two- 
dimensional dynamical system, and according to a the- 
orem by Poincare |f lj it cannot exhibit such behaviour. 
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FIG. 1: The direction of the J-flow given by Eq. (|5| in the 
bistable regime. Also shown are two sample integral curves, 
one (in green) with initial conditions causing it to spiral into 
the low-amplitude attractor, and one (red) which ends up 
in the high-amplitude stationary state. The third irregular- 
ity in the flow, which can be discerned near (0.500,-0.496), 
corresponds to the saddle point. Integrating the flow back- 
wards in two directions starting out from the saddle, one ob- 
tains the separatrix demarcating the two basins of attraction, 
which spirals out to infinity in an appealing Yin- Yang pat- 
tern, shown in Fig. |10| Physically, while all initial states near 
the origin are attracted to the low-amplitude state, the fate 
of initial states with large initial amplitude depends on their 
phase, so that the two basins of attraction continue to be 
interleaved for arbitrarily large radii. 



The difference is of course the specialization to infinites- 
imal e in the latter. 

The stationary states, which were limit cycles [TT] of 
the original Eq. Q, become fixed points in the rotating 
frame upon coarse-graining, and can be found by setting 
the right-hand side (rhs) of Eq. Q to zero. In polar 
coordinates u — A cos <f>, v — A sin 0, one obtains }13) 



cot ( 



' 2 +Qa'A' 2 -Q') 2 A' 2 
a' A' 2 - Of , 



(7) 
(8) 



where A' = AcJoj/f is defined similarly to U,V in 
Eq. (J6J). For a' i 0, these reduce to the Lorentzian re- 
sponse of the damped harmonic oscillator, but even in 
the nonlinear case, the driving / merely sets a scale for 
the amplitude A. As a result, the amplitude response 
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Rescaled Amplitude vs. Rescaled Detuning at different f 




FIG. 2: The reduced amplitude response of the Duffing os- 
cillator. The low- (red) and high-amplitude (blue) stationary 
states, as well as the unstable branch (green), are the classical 
response Q. Data points denote the stationary quantum re- 
sponse at T — 0, see below Eq. (21 1; the stronger the system 



is driven, the closer it is to the classical limit, and the steeper 
the curves are near the crossover point Q' w 3.3. 



A'(Q') is a one-parameter family of curves, dependent 
only on a' . For a' > (a' < 0) these curves become 
asymmetric and skewed to the right (left), until at the 
critical value a 1 = 32/3 5/2 « 2.0528 the solutions bifur- 
cate, and become multiple- valued for an increasing range 
of detunings Q' — the bistable window (tri-colored line in 
Fig. |2|. On the one hand, a choice for a' should ex- 
ceed the critical value sufficiently to afford a substantial 
bistable window. On the other, if a 1 is taken too large, 
the separation between the two stationary states becomes 
so large that the quantum tunneling phenomenon dis- 
cussed below is too weak to be observed. For ease of 
comparison, in all our examples we take of — 6, with 
bistability for 2.8568 < ft' < 4.5563. Finally, note that 
the maximum amplitude is always A' = 1, attained for 
il' — |a', as evidenced in Fig. [2] 

Quantum theory. Thus far, the Duffing bistability is 
analogous to the one of a particle in a double- well poten- 
tial, say 



U(x) — x 4 — bx 2 + cx 



(9) 



For c = 0, the unique minimum U(x=0) = for negative 
b bifurcates at the critical value b — 0. For b > 0, there is 
an increasing (and in this case, symmetric) range of bi- 
ases c for which U{x) features two minima, until the left 
(right) minimum loses local stability at the lower (upper) 
end of the bistable c-window. Thus, the parameters b and 
c play the roles of a' and f2', respectively. Quantum me- 
chanics in a potential such as ^ has several outspoken 
features. For instance, for b > the ground state rapidly 



changes character as a function of c in a narrow range 
(of the order of the tunneling amplitude, hence exponen- 
tially small near the classical limit) around c = 0, from 
being localized in the left well, via being a symmetric su- 
perposition, to being localized in the right well. We call 
this phenomenon crossover; for U(x) in (19]), the crossover 
point c = is obvious by symmetry, but even for more 
complicated U, one merely has to find the bias point for 
which the well minima have the same potential. 

Guided by the above analogy, in the quantum case one 
expects tunneling transitions between the two stationary 
states to become possible also for the Duffing model. Key 
questions are to which extent such transitions resemble 
the (equilibrium) tunneling in a potential such as (|9| on 
the one hand, and the classical dynamics of Eqs. ([T])-([8]) 
on the other, in particular close to the classical limit. 
Crossover phenomena are of interest as well, especially 
since the location of the crossover point does not seem to 
be a priori obvious. In Ref. 21 the focus is on the JBA 
application of a quantum Duffing oscillator coupled to a 
qubit. In Ref s. [S] and® quantum Duffing dynamics was 
studied near the edge of the bistable window (called the 
critical case in the former and the bifurcation point in 
the latter) . In Ref. \7\ the quantum Duffing oscillator is 
studied using a superoperator Floquet technique. 

Quantum tunneling in driven systems has been stud- 
ied in many contexts |14j . but here the very location of 
the stationary states depends on the damping in leading 
order (since A ~ r y^ 1 A', with A' of order one). As a 
result, the problem does not have a weak-damping limit 



[indeed, the evolution equation (111 below has damping 



coefficients equal to unity in its last two terms] , and a de- 
scription in terms of the density matrix p must be used 
from the outset. 

We use a quantum Master equation (QME) formula- 
tion. The use of a (Markovian) QME for damped quan- 
tum oscillators is not always straightforward, or indeed 
uncontroversial ISj. However, it is justified here since 
(a) the damping ~ is very weak compared to the har- 
monic motion, and (b) we consider only the long-time, 
coarse-grained dynamics. Note further that since also 
the nonlinearity is ~ e, any cross-effect between it and 
the dissipation is at least 0(e 2 ), and beyond our consid- 
eration. This justifies the use of the standard harmonic- 
oscillator damping terms [TB] also for the Duffing oscil- 
lator. Because of the restriction to long times, we work 
exclusively in the rotating frame 



p{t) = U\t)p{t)U{t), U(t) = 



jNt 



N = a^c 



(10) 



which are analogous to the classical Eq. Q. Upon per- 
forming the coarse-graining or rotating-wave approxima- 
tion (RWA), one obtains without further ado the QME 
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in Lindblad [T7] form 

i— — = [H, p] + -(n+l)(2apa i — a)ap — pa 1 a) 

UT 2 

+ -n(2a^ pa — aa) p — paa^) 
= C{p} , 



(11) 



(12) 



where we introduced the Liouville superoperator C and 
the scaled rotating-frame Hamiltonian 

H = -^N+£=(a + a^ + l^(N 2 + N). (13) 

Here, the last term oc a' denotes the potential nonlin- 
earity; the form iV 2 + N arises by rearranging creation 
and annihilation operators in the expansion of the po- 
tential x 4 , retaining only co-rotating terms [cf. above 
Eq. Q]. The time scale r is as in Eq. Q. Compared to 
the classical case, however, there is an extra parameter: 
the dimensionless driving 



/' = 



7 V 



(14) 



Namely, while the classical energy can be scaled away, 
in the quantum case the photon number TV is discrete, 
with N 3> 1 representing a semiclassical situation while 
N < 1 is the deep quantum limit. This is reflected in 



the definition (14), from which one sees that the classi- 
cal limit /' — > oo can be achieved either through strong 
driving /, taking a large mass m, or letting h — > 0. We 
have also introduced the Bose occupation number 



1 



(15) 



so that for T — > 0, only the first (decay-type, ~ a pi) 



dissipative term survives in Eq. ( 11 ), while the excitation 
term ~ cl* p vanishes. 

For initial exploration, and confirmation of quantum 
tunneling, Eq. (11) can be simulated directly by trun- 
cating N, choosing a p(r=0), and using a standard RK4 
solver for the resulting ODE system. Care has been taken 
to ensure probability conservation for arbitrary trunca- 
tion and discretization of the time steps; however, the 
exponentially small tunneling effects are very sensitive 
to cutoff artifacts, and convergence must be verified sys- 
tematically. A representative example is shown in Fig. |3j 
starting from the vacuum state p~(0) = |0)(0|. Initially, 
the driving force rapidly increases the photon popula- 
tion, until it becomes counterbalanced by the dissipation 
(which is more effective for higher states). For times 
t ~ 10 (top panel), the system seems to settle in a wave- 
packet type state, which is the quantum counterpart to 
the classical low-amplitude state (see Fig.|2]for quantita- 
tive verification). However, if one persists in simulating 
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FIG. 3: Real-time evolution in the rotating frame on the long 
(top) and ultralong (bottom) time scale, showing the proba- 
bility distribution in the number basis as a function of time, 
for detuning Q' = 3.24 and driving /' = 7. While only the di- 
agonal elements of {n\p(r)\m) are thus shown, of course these 
data have been obtained from a simulation of the full density 
matrix. 



until r ~ 10 3 (bottom panel), an entirely new feature 
emerges in the form of a wave packet corresponding to 
the classical high-amplitude state, which becomes popu- 
lated by quantum tunneling from low amplitudes. Only 
for r ~ 10 4 is a genuine steady state reached (we prefer 
to avoid the term "equilibration" for this driven system) , 
through a balance between low-to-high-amplitude tun- 
neling and vice versa. Thus, nonequilibrium tunneling 
has opened up a new, ultra-long time scale in the quan- 
tum Duffing model. 

Given that the RWA has yielded a time-independent 
evolution operator, simulation is ultimately less power- 
ful than directly extracting the r —¥ oo tunneling dy- 
namics by studying the eigenvalues of C in Eq. (12) [7]. 
We have performed this by numerical diagonalization, 
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using a non-hermitian sparse- matrix routine [T5]. The 
main disadvantage of this method is size: while p is 
N x N, the matrix form of C is N 2 x TV 2 . With some 
optimization, we have used truncations up to N ~ 120 
and substantially larger N would have been possible, if 
very time-consuming. While this method thus is slightly 
laborious, especially close to the classical limit where 
quantum-classical correspondence must be studied, we 
have found its excellent stability and reliable convergence 
to be worth the effort. 

As a result, one obtains eigenvalues Xj and eigen- 
vectors" pj, in terms of which the time evolution can 
be expressed as 



P(r) 



(16) 



In contrast with Hamiltonian spectra, only the sum p(r) 
is required to be physical, while the individual fij need 
not have all the properties of density matrices, and can 
for instance be non-hermitian. From the structure of C 



in Eq. (11), one readily deduces that for each eigenpair 



(\j,Pj), one also has the pair (— A*,/^), and this sym- 
metry is indeed observed in the numerical output. Thus, 
the eigenvalues either occur in conjugate pairs, or are 
purely imaginary. Further, it is essential for stability 
that no eigenvalues lie in the upper-half complex plane, 
which would lead to a divergent evolution oc e 



Im XjT 



Eq. (16) 



We sort the output of the iterative diagonalization by 
increasing — Im Xj ; clearly, the more of these one wishes 
to represent accurately, the more photon-number states 
must be included. With this convention, we find that in 
fact ImXj < for all j, with the exception of 



Ai = 



(17) 



representing the unique stationary state. Further, with 
increasing /', one observes a separation of time scales il- 
lustrated in Fig. [4] while Aj>3 remain of order one, A2 
(which thus is purely imaginary) becomes exponentially 
small, so that it must be carefully resolved from Ai. Phys- 
ically, on the time scale r ~ 1, the quantum system re- 
laxes towards quantum counterparts of the two classical 
states within each basin of attraction, analogous to the 
dynamics shown in Fig. [I] On the new, ultra-long time 
scale ~ I A2 1 1 , quantum tunneling exchanges population 
between these two states, until the true stationary state 
Pi is reached for r — > 00. As a consequence, for r> 1 
one has, with exponential accuracy, 



p(r) « p x +c 2 p 2 e" |A217 



where only c 2 depends on the initial conditions, 
bility conservation implies 



Trp 2 = , 



(18) 



Proba- 



(19) 



E 




FIG. 4: In the bistable region, only the eigenvalue A2, which 
is related to tunneling, decays exponentially with the driving 
strength /'. This behaviour is more pronounced near the 
crossover point. 
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FIG. 5: Diagonal elements of the stationary density matrix 
pi (a) and the tunneling "density matrix" p*2 (b), for detuning 
f2' = 3.2 4 and driving /' = 9. In the two-state approxima- 
tion ( 20 1 , p2 contains no new information, since its peaks are 



the same shape as those of pi, but are constrained to have 
equal area and opposite sign (the overall prefactor of p~2 is 
arbitrary). Instead, the dynamical information is contained 
in the eigenvalue A2, see Fig. [4] 



which is indeed satisfied by the numerical result, see 
Fig. [5} In practice, Eqs. (17) and (19) serve as a useful 



test for convergence with respect to the number of in- 
cluded photon states (the details of this convergence are 
sensitive to the exact manner in which C is truncated). 

The magnitudes of {n\p\\m) are shown in Fig. [6j for 
parameters in the bistable regime close to the crossover 
point. The two lobes centered on (L,L) and (H,H) are 
readily distinguished, with L w 6 and H ss 33. Their ex- 
tent away from the main diagonal shows that p\ contains 
two semi-classical states with significant partial coher- 
ence (the quantification of which is possible, but not our 
priority here). We have examined these same data also 
on a logarithmic scale, down to machine precision. While 
the lobes of Fig. [6] are then seen to have tails extending 
far away from the main diagonal, no further structure 
is observed centered on (L,H) and (H,L). Thus, while 
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FIG. 6: The magnitudes of the stationary density-matrix ele- 
ments in the number basis (n\p\\m), for f2' = 3.23 and /' = 9. 



each lobe is partially coherent, p\ is a mixture of them 
and not a superposition. Together with A2 being purely 
imaginary, this supports our claim that the type of tun- 
neling discussed here is incoherent — not surprising given 
the prominent role of dissipation in the Duffing oscilla- 
tor [H]. 

The above implies that, for the longest times, the relax- 



ation ( 18 ) is governed by a simple two-state rate equation 



± (Ph 
dr \P L 



fPn 
\Pl 



(20) 



where denotes the transition rate from the high- 
amplitude semiclassical state to the low-amplitude one, 
etc. Equation ^ implies P H = 1 - P L = r T /(r r +r i ) 
for the stationary distribution and IA2I = + T^ for 
the relaxation rate. It follows that A2 is dominated by 
the larger of the two transition rates, so that it is ex- 
pected to be small only in the middle of the bistable 
window. Namely, near the boundary, the upper well in 
the potential-tunneling analogue ^ becomes shallow, so 
that the corresponding escape rate is large; see Fig. [4] 
Thus, IA2I should attain its minimum for, or at least close 
to, r-j- = T±, in which case Ph = Pl = \ — the crossover 
point. Properly, the model (20) should be regarded as 



phenomenological, and pi and | A2 1 are independent nu- 
merical data. Either equiprobability for the former or 
the minimum in the latter can be used to pinpoint the 
crossover, with slightly different but consistent results 
(black lines in Fig. [8]). 

The correspondence to classical Duffing dynamics is 
vividly illustrated by the amplitude response (data points 



in Fig. [2]). To determine the proper expression for A' , one 
can transform p\ back to the lab frame using the inversion 



of Eq. (10) and calculate (x}(t) [the t-dependence arising 
solely from the one in U(t)]. An unexpected feature is 
the small "bump" visible on the low-amplitude side of the 
crossover; like the remainder of the crossover curves, this 
bump becomes narrower closer to the classical limit, but 
it retains a finite amplitude and never disappears. The 
issue is clarified by noting that in our driven system, the 
phase angle with respect to the driving force is relevant, 
and therefore one should properly introduce a complex 
amplitude Z = Ae 1 ^ . Using the same considerations as 
above, one finds in dimensionlcss variables: 



Z' = I £ ^W)(n+l\~ Pl \n) = ^ Tr[p ia ] 



r 



(21) 



indeed, A' in Fig.|2]was simply calculated as (A') = \Z'\. 
Thus, the amplitude crossover properly takes place in the 
complex plane, see Fig. [7| Amusingly, following the ap- 
proximately straight crossover line formed by the data 
points from the low- to the high-amplitude branch ini- 
tially brings one closer to the origin, explaining the non- 
monotonous behaviour of (^4') (O'). Indeed, in a crossover 
plot of, say, the energy ex (N), no such bump is observed. 
To be sure, this "complex phase" phenomenon has noth- 
ing to do with quantum mechanics, and the amplitude 



response in the case of classical diffusion ( 22 ) below fea- 
tures the "bump" as well. 

The approach to the classical response, where appro- 
priate, provides one means of verifying our numerics. Di- 
rect comparison with analytical results is limited to the 
harmonic case oi — 0, in which C can be diagonalized us- 
ing superoperator Lie algebraic techniques [20] . We omit 
the details, since this case is more easily dealt with in the 



coherent-state formulation ( 33 )— ( 34 ) below. However, we 



do note that with further specialization to zero temper- 
ature, the ensuing stationary state p\ is pure even in the 
presence of dissipation, namely, a coherent state. The 
reason is that coherent states are eigenstates of the anni- 
hilation operator a, which generates the only nonunitary 
process at T = 0. All these are in full agreement with 
our numerical results. 

Classical diffusion at finite temperature. The ultra- 
long tunneling time scale discussed above is only "with- 
out classical counterpart" at temperature T = 0. At 
T > 0, noise-assisted transitions from one stationary 
state to the other are expected to occur, in analogy to 
thermal activation across the barrier of bistable poten- 
tials such as ([9]) [21 . Thermal and quantum processes 
can be studied in combination by analyzing Eq. (11) for 



finite n; for now, we focus on thermal effects in the purely 
classical case [22]. 

Extending the Newton-Duffing dynamics into a dif- 
fusion process is most reliably done in the lab frame, 
by adding a Langevin term ^/le^k^T '/m£(t) to the 
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FIG. 7: Amplitude crossover in the (U, V) plane. Classical 
and quantum (for /' = 9) data are parametric plots for 2 < 



Q' < 5, the former using Eqs. |7|) and 
using pTl for Z' = U + iV [cf. above Eq. 



and the latter 
1 . Note that all 



classical stationary states fall on the circle TT^+ (V+ 1) 
as is readily proven. 



l 

4 ' 



rhs of Eq. ([IJ, with £(t) being (5-correlated Gaussian 
white noise. Equivalently, one can make the fluctuation- 
dissipation substitution ejd p p e^d p {p + mk^Tdp) in 
the associated Liouville (classical phase-space continuity) 
equation [23] . 

Transformation to the rotating frame can be done in 
either formulation, with equivalent results. In Langevin 
language, the p-noise in the lab frame becomes a rapidly 
rotating noise term in Van der Pol coordinates, which 
should be coarse-grained like the other evolution terms, 
keeping in mind that the noise autocorrelation time is the 
only time scale which is short compared to our "fast" 
scale ~ uj^ 1 . This results in separate, uncorrelated u- 
and v- noises to be added to Eq. Q. Alternatively, the 
diffusion term oc d 2 can be transformed using Eq. |3j), 
upon which the coefficients of the resulting diffusion ma- 
trix are time-averaged separately. In either case, one 
obtains the Duffing-Fokker-Planck equation (FPE) with 
isotropic diffusion term for the classical probability den- 
sity p(U,V), 



= 4-p{p} . 



(22) 
(23) 



where V = (d(j,dv) T denotes the gradient with respect 
to the scaled coordinates of Eq. ([6]), and with «7, r as 
in Eq. Q. Instead of the reduced temperature feeT/fi^o 



entering Eq. (11) through Eq. (15), here one encounters 
the "classical temperature" T c \ = ksTj 2 /mf 2 , compar- 
ing to the energy 
with amplitude A' = 1. 



muj 2 A 2 of classical oscillations 



Equation ( 22 ) can be studied by expanding it in a dou- 



ble Hcrmite basis [2"4"j . 

p(U,V) = J2 C nm€ ) m^(V). (24) 
nm 

Even though physically there are no creation and annihi- 
lation operators in classical diffusion, quantities such as 
U, djj , etc are readily represented in matrix form using 
the recursion relations for Hermite functions, so that C„„ 



in Eq. (23) is sparse similarly to C in (12). The result 



ing numerics are strongly analogous to those resulting 
from Eqs. ( [TT| ) (13) in the quantum case, and require 
resources of the same order [25] . As T c \ is lowered, the 
distribution p(U, V) retreats to two narrow peaks cen- 
tered on the stationary states (colour-scale component 
of Fig. 10), requiring increasingly many basis functions 



to represent accurately. Normalization of p involves the 
doubly even coefficients C2 n ,2m 2(i . Extraction of the 
mean amplitude (A') = (\/U 2 + V 2 ) and of the prob- 
ability distribution over the two peaks (for sufficiently 
low T c i that such a quantity is well-defined) is compar- 
atively less elegant, and can be performed on a grid by 
transforming back from matrix form to real space using 



Eq. (24) 



Thus, the classical finite-temperature amplitude re- 
sponse and relaxation rate are obtained. As a function 
of decreasing T c i, these are found to behave qualitatively 
similarly to the T — quantum data of Fig. [2] as a func- 
tion of increasing /'. However, comparing the two sets of 
data in Fig. [8] reveals a curious difference: while each set 
individually is internally consistent (in that the ambigu- 
ity in the definition of the crossover point vanishes as the 
width of the crossover tends to zero) , there is a small but 
significant difference in the value of the crossover point 
predicted by each set. 

Since the only difference is that the QME data are at 
T = for h -> while the FPE data are at h = for 
T — > 0, the system's behaviour in the h — > 0, T — > limit 
seems to depend on the order in which these two lim- 
its are taken |27j . If confirmed, this represents a strik- 
ing difference from equilibrium tunneling in potentials 
such as Q, where crossover always occurs at potential- 
minimum degeneracy (even for asymmetric potentials, 
since zero-point motion in the well minima etc all van- 
ish in the low-T classical limit). The remainder of this 
paper is devoted to independent verification of this phe- 
nomenon, as well as to understanding of its origin. 

Semi- analytical low-T asymptotics. For an indepen- 
dent verification of what was obtained above through ex- 
trapolation of large-scale numerics, we have performed 
an asymptotic analysis. For the FPE (22), this means 



attempting to take the low-T limit analytically, through 
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FIG. 9: Equation (26 1 constrains it to a semicircle with di- 
ameter —J. 
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FIG. 8: The crossover point of the Duffing model, as extracted 
from finite- T classical (red) and T — quantum data (black). 
In either case, the dashed line links data points obtained us- 
ing a (static) equiprobability condition, while the solid lines 
are obtained if one uses the minimum in the (dynamic) relax- 
ation rate [A2I to pinpoint the crossover. Extrapolating any 
curve to the j/-axis realizes one manner of taking the h, T — > 
limit. For each color, dashed and solid lines extrapolate to the 
same Q! in the limit of abrupt crossover. However, the red 
and black curves do not agree with each other. The red star 
on the t/-axis refers to the classical crossover point found in 
Fig. 11 using data obtained from Eq. (|28[l, in which the T — > 



extrapolation is performed analytically. 

a quasi-Boltzmann Ansatz 

p(U,V,T cl ) ~ e-^ u < v ^ , 



(25) 



which plays a role here analogous to semi-classical op- 
tics or the WKB approximation in their respective fields. 



Substituting into the stationary limit of Eq. (22 1, to lead- 
1 one obtains 



ing order in T cl 

R (R+ J) = 



R=VS 



(26) 



Thus, the original linear second-order PDE (22) has been 
reduced to a nonlinear first-order one. 

For a gradient flow J 
solved by S — U, 



this leads to 



[2R{J,x) + J]-V X 



dJu 


dJv ) 


dV 


dU ) 


/dJ v 


dJ v 


\ dU 


+ 

dV 



(1 - cosx) 

(28) 



sm-X 



Through the method of characteristic curves [25], this 
PDE is equivalent to a system of three coupled ordinary 
differential equations (for U, V, and x), to which one 
may add the integration of VS 1 = R( J, x) as a fourth for 
simultaneous solution. 

The resulting differential system is highly nonlinear 
and unstable, with a parasitic solution x = =>■ R = 0. 
However, it has the key feature that the low-T limit is 
already incorporated. To integrate a solution for S one 
needs an initial value, which can be obtained at the sta- 
tionary points of J, near which to leading order 



Ju 
Jv 



a b 
c d 



dU 

sv 



(29) 



with different coefficients a through d for each of the three 
stationary points. Since one can convince oneself that 
these points are stationary also for R and for \, the rhs 



of (28) can be set to zero there, yielding 

+ d 



tan 



(30) 



-VZ7, Eq. (26) is immediately 
but for J as in Eq. (5L the analy- 
sis is more involved. First, we solve Eq. (26) point-wise 
through the parametrization 

R= ~ + ^-cosx+ ^sinx , (27) 

with J± = (—Jv,Ju) T (see Fig. [9J. Instead of two un- 
known components of R, one then deals with a single 
scalar field x(U, V). To ensure that R can be uniquely 
integrated into an effective potential S(U, V), one has to 
impose Vxfl=0 explicitly. Taking the curl of Eq. ( 26 1 , closely; a representative case is shown in Fig. 10 



One cannot start at a stationary point since there would 
be no evolution, but one can start close to it, taking 
(5U,SV) = (cos0,sin0X with C ~ HT 9 -l(r 6 (check- 
ing for convergence with respect to further decreasing Q . 
The ultimate aim is to compare the values of S in the 
low- and high-amplitude "wells" , which is achieved by a 
shooting method, tuning the angles 9 such that one ob- 
tains curves linking either well to the saddle point, where 
for definiteness we set S = to fix the immaterial overall 
additive constant in S. Implementing the method, we 
have found that different trajectories shot from the sad- 
dle to the wells converge exponentially, leading to loss 
of resolution. Thus, we start from the wells, and shoot 
out x-curves for a range of values for 6. Since the saddle 
is an unstable stationary point of these curves, must 
be tuned extremely accurately to approach the saddle 
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FIG. 10: Shooting to the saddle point from the vicinity 
of the low- and high-amplitude stationary states, at detun- 
ing Q' = 3.3. Red/yellow: classical particle trajectories, in- 
tegrated from the flow in Fig. [T] Black: separatrix. Ma- 
genta/brown: x-curves from Eq. ( |28[ ), aimed towards the sad- 
dle to determine the effective barrier height. Background and 
colour scale: effective potential approximated via — In pi [cf. 
Eq. |25f] for T cl = 0.0035. 



According to the above analysis, if the ^-curves are 
indeed part of a global solution S(U,V), then Eq. (30) 



should again be satisfied at the saddle (which has its 
own coefficients a through d). We have confirmed that 
this is indeed the case, for both low- and high-amplitude 
starting points; however, Eq. p0| ) is not satisfied at other 
points along the curve. Technically, the v-curves play 
a role similar to instanton paths in the WKB method; 
however, in the present setting we do not have a physical 
interpretation of these curves. 

Carrying out the shooting procedure for a range of Q' 
in the crossover region, we acquire the data shown in 
Fig. [IT] Reconstructing the effective potential S(U, V) 
gives us our closest analogy with the potential bistability 
of Eq. ((9|. The inferred crossover point f2' s» 3.292 is 
shown as the red star in Fig. [8] The excellent agreement 
with the extrapolated diagonalization data is compelling 
evidence that the latter indeed have the claimed accuracy 
(at least for the classical case). Note also that Fig. [IT] 
predicts the high-amplitude state to dominate at lower 
detuning and vice versa, as is indeed observed in Fig. [2j 
Not only the crossover point, but also the correspond- 
ing barrier height AS is instructive. Namely, from an 
activation-energy consideration of Eq. ( 25 ) , one expects 
— In |Aa| ~ 2AS/T C \ + const. Thus, AS* can alternatively 
be inferred from the relaxation data \2(T C \) in a most 
traditional analysis the Arrhcnius plot 29J. We have 
verified that this yields AS" w 0.018 ± 0.002, in agree- 
ment with Fig. 
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Thus, static (pi) and dynamic (A2) 
data lead to consistent conclusions, as was also the case 
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FIG. 11: The effective potential S versus detuning for the 
two classical stationary points. The two curves intersect at 
the crossover point. 



for the QME when inferring the incoherent nature of the 
tunneling either from Re A2 = 0, or from the absence of 
off-diagonal peaks in Fig. [6j 

Coherent- state analysis. The above successful asymp- 
totic study of the FPE p2l prompts one to attempt 



a similar (large-/') analysis for the QME (11) as well. 
However, the number (Fock) basis is notoriously un- 
suited for studying semiclassical properties of oscillator 
systems. Instead, it is logical to employ the "most clas- 
sical" oscillator states — the coherent states (in fact sim- 
ply Gaussian wavepackets in the position basis) \Q = 
e~ K|2/2 IT=o y=rN' with ( n > = Kl 2 - Due to their over- 
completeness |(H0| 2 — e"'""^ 2 , one has not one but two 
coherent-state representations: 



P = /d 2 CP(C)|CXCI 



Q(C) = -(CIpI0 

7T 



(31) 
(32) 



Either one has its strengths and drawbacks. From 
the above definitions, one can immediately calculate 
P i->- p n- Q. For our diagonalization data, Q(U, V) 
is shown in Fig. |12| since it is real, such a plot in 
principle characterizes the system's state fully, unlike 
Fig. [6j where the phase of the (off-diagonal) density- 
matrix elements is suppressed. Indeed, comparison of 



Figs. [10| and 12 makes the classical-quantum correspon- 
dence manifest. The reverse operations Q > p > P are 
more subtle. Some plausible-looking Q are simply un- 
physical (i.e., not corresponding to any p), for instance 



since Eq. (32) endows Q(C) with a minimum linewidth. 
Also, while a P-function can be proven to (uniquely) ex- 
ist for any p |16j . such functions are highly singular even 
for some simple p's (e.g., pure number states). 

Almost by definition, a\() — £|C), and one also veri- 
fies at |0 = (<9 C + C*/2)|C). Thus, Eq. |uj| can be tran- 
scribed term by term into coherent-state language. In 
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FIG. 12: The coherent-state Q-function, obtained from its 



definition ( 32 I for the stationary state pi , close to the (T = 0) 
crossover point. 



carthesian coordinates (Re£,Im£) = (f'/y/2)(U, V), the 
results read 



dP 



dQ 
dr 



= -v • (JP) 



3a' 



If" 



r V 2 P 



+ ^2 [i d u - 9 2 y)UV - d v d v {U 2 - V 2 ) 

+ 2(d v U -d v V)]P 
-V • (JQ) 



(33) 



n+1 9 „ 
^V 2 Q 



V' 1 



3a' 
8f 



^[{d 2 ~ d 2 v )uv - d v d v {u 2 -v 2 ) 

+ 2(d v U - duV)]Q . 



(34) 



The "quantum terms" take a transparent form in polar 
coordinates, 

(dlr - 8l)UV - d v dv(U 2 - V 2 ) + 2(d v U - d v V) 

= -d A ,d^A! . (35) 



From the definitions (31 1, (32) it follows that the 



Q-function is merely a Gaussian smoothing of P, viz., 
Q(() = I fd 2 v P(v)e-\ u ~ c \ 2 . Through integration by 
parts, this can be used to check the consistency of Eqs. 



(33) and d34 



In the harmonic case a' = (note that a' also occurs 
implicitly in J), Eqs. (33) and p4| ) have a readily found 
Gaussian solution 



P(U, V) = — exp 
nn 



r 

2n 



n(u,v) 



H(U,V) = \ u- 



1 + 0' 2 



V 



1 + 0' 



(36) 
(37) 



with n n- n + 1 for Q(U, V). In particular, for n — > the 



P-function tends to a 5-peak, which according to Eq. (31 ) 



corresponds to a coherent state [cf. below Eq. (21)] 



The coherent-state equations ( 33 ) , ( 34 ) do not readily 



translate into a numerical scheme. While they would 
be satisfied by numerical solutions for P and Q once 
found by other means, forward integration in r is un- 



stable since the "quantum terms" ( 35 1 correspond to an 
indefinite "diffusion matrix" , as opposed to the positive 
definite V 2 . More properly, the rhs of Eq. (35) demon- 



strates that these terms represent a hyperbolic operator, 
as opposed to the elliptic V 2 [2"5] . 

Instead, the strength of these equations is in making 
the various limits of interest more transparent. First of 
all, consider the classical limit at finite (classical) tem- 



perature. Since n/f « T c \ for n >• 1, both Eqs. (33) 



and ( 34 ) converge to Eq. (22 1 if /' — > oo at fixed finite T c i, 
the difference between them becoming negligible. Thus, 
classical diffusion is almost trivially obtained as a limit 
of the dissipative quantum dynamics. We have verified 
this behaviour numerically for both the stationary ampli- 
tude response and the tunneling rate (data not shown), 
using Eq. (15) to determine the n(f',T c {) to be used in 

Eq. ^TTJ> 

Contrast the above with the limit /' — > oo at fixed n 
(implying T c \ j 0). In this case, both the (diffusive) 
V 2 -term and the "quantum terms" are 0(f'~ ), so that 
neither dominates the other and both must be retained. 
The most striking case of this type of limit is n = 0; for 
the P-function, the diffusive term then is absent entirely, 
and only the "quantum terms" cause a deviation from 
the classical contractive flow of Fig.[l] Thus, without ac- 
tually having solved them, Eqs. ( [33] ) and (34) do yield the 
definite prediction that (and why) the limit /' 1 , T c \ — > 
depends on the value of n, as was observed by various 
means above. 

In an attempt to quantify these findings, for n = 
we have pursued the asymptotic Ansatz analogous to 



Eq. (25) 



P(U, V, /') ~ exp 



3a 



r H(U,V) 



(38) 



Since it was pointed out above that the second-order dif- 
ferential operators in Eqs. (22) and (33) have an essen- 



tially different character, and since we do not have inde- 
pendent numerical data for P, such an approach is admit- 



tedly heuristic. Substitution of Eq. (38) into (33) leads 



to a considerably more complicated analysis than the one 
leading up to Eq. (28) in the diffusive case — partly be- 



cause of the position dependence in Eq. (35), and partly 



because the counterpart to Eq. (26) constrains VH to 
hyperbolae not circles. We omit the details, since these 
considerations have been superseded by an unexpected 
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FIG. 13: Crossover point for n — (i.e., low-temperature 
limit first) as determined by degeneracy of the "Quantum 
potential" p9i. 



analytic solution: 



H(U,V) 



u 



u 2 + v 2 



with <f> the phase angle as above Eq. ([7|. 

One verifies that Eqs. (38 1, (|39|) solve the n 



(39) 







coherent-state equation (33) not only to leading order 



in f~ , but in fact exactly. Construction of a crossover 
plot analogous to Fig. [TT] merely involves the evaluation 
of the "quantum potential" T-L at the classical stationary 
points, again using its value at the saddle point to fix 
the additive constant. The resulting data are shown in 
Fig. [131 

It must be repeated that the whole procedure start- 



ing with Eq. (38) is heuristic. We have not proven this 



asymptotic form for P; given that the ensuing H has at 
least two major unphysical features (its multiple-valued 
character due to the occurrence of <f> and its singularity at 
the origin), the Ansatz may indeed not be wholly correct. 
Nonetheless, the agreement of its upshot O' ps 3.222 with 
the extrapolated numerical value for the n = crossover 
point is unlikely to be mere coincidence. For the mo- 
ment, we have to leave the solution's precise status as a 
subject of further study. 

Conclusion. Let us review what has been accom- 
plished. We have the two states visible in Fig. [I] whose 
definition requires nothing beyond 19th-century physics. 
The question which state is "more stable" (or equiva- 
lently, for which parameter these states are "equally sta- 
ble") could be expected to be answerable in a similar 
setting. However, this intuition has been confounded; 
instead, the answer essentially depends on the details of 
how this limiting classical picture is reached from the un- 
derlying finite-temperature quantum theory. In contrast, 



for equilibrium bistability described by Eq. ([9]), the corre- 
sponding determination can be made immediately from 
a plot, keeping in mind the word's literal meaning "scales 
in balance" . 

It should be investigated whether other instances of 
nonequilibrium bistability |31j exhibit analogous phe- 
nomena and if so, whether a general description is possi- 
ble. In any case, similar subtleties occur in several places 
in quantum physics; see, e.g., Ref. [32] for a case of non- 
commuting zero-damping and low-temperature limits in 
the Casimir effect. Given that one needs extremely long 
times to reach the stationary state in the crossover regime 
[~ (eA2) _1 , with A2 as in Fig. [4], one could verify whether 
going beyond the RWA affects any of our conclusions. 

Moving beyond the specifics of the nonunique 
crossover, our investigations have yielded extensive data 
and calculation schemes for quantum Duffing oscillator, 
which is relevant to the qubit-readout application [3J. 
Acknowledging previous efforts [IHZj, the present results 
may be the most accurate and general [in the chosen 
regime of small e in Eq. Q]. 
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